% Copyright 2009 Neurosciences Research Foundation, Incorporated
function [cost_measure spikeout mfr all_target_mfr] = spnet( conn_data_in, animate, optional_pattern, all_target_mfr)

global s sd syntype D N syndata 
global  Ne IN A AI
global r vhist uhist ehist ampahist nmdahist gabaahist gababhist itothist se Ie

N_connections = 6;
conn_vars = 1;
N_groups = 3;
group_vars = 0;


draw = false;

%search_params = [1 3 4];
search_params = [3];
conn_data = zeros(N_connections,18);


if nargin < 2
    draw = false;
    animate=false;
end



noise = 0.2;
signal = 1.0;


% excitatory neurons   % inhibitory neurons      % total number 
Ne=100;                Ni=50;                   N=Ne+Ni;

% set up group sizes.
IN = 1:Ne/2;
A = max(IN)+1:Ne;
AI = max(A)+1:N;
%N = max(AI);


conn_name{1}='IN->A';
conn_name{2}='A->A';
conn_name{3}='AI->A';
conn_name{4}='IN->AI';
conn_name{5}='A->AI';
conn_name{6}='AI->AI';


if exist('conn_data_in') && ~isempty(conn_data_in)
    if size(conn_data_in,2)==1
        %group_data_in = reshape(conn_data_in(N_connections*conn_vars+1:end), N_groups, group_vars);
        conn_data_in = reshape(conn_data_in(1:N_connections*conn_vars), N_connections, conn_vars);
        if draw
            group_name{1}='IN';
            group_name{2}='A';
            group_name{3}='AI';
            
            
%             disp('          ');
%             for i=1:N_groups
%                 fprintf('%10s ',group_name{i});
%                 fprintf('%10.3f ',group_data_in(i,:))
%                 fprintf('\n');
%             end
            disp('          std        wt range,      radius  background ampa, gaba,   nmdagain, gababgain,      f,     d,       u')
            for i=1:N_connections
                fprintf('%10s ',conn_name{i});
                fprintf('%10.3f ',conn_data_in(i,:))
                fprintf('\n');
            end
        end
    end
    
else
                      % background ampa, gaba,   nmdagain, gababgain,   f   d u
                      % (first 2 are background spike mfr for ampa inputs and gaba
                      % inputs.}
%     group_data_in(1,:)=    [0 0 .5 .1 5 300 .5];
%     group_data_in(2,:)=    [800 400 0.5 .1 1000 800 .5];
%     group_data_in(3,:)=    [0 0 .5 0 20 700 .2];
end
start_save_sec = 0;
stop_sec = 1;



gampa_background=zeros(N,1);
ggabaa_background=zeros(N,1);



epoch_ms = 1000;
npatterns = 4;
isi_ms = floor(epoch_ms/npatterns);

if nargin > 2
    r=optional_pattern;
    
    w=cos(-1.5*pi:1.5*pi/(Ne/4):1.5*pi);
    w(w<0)=0;
    for i = 1:npatterns
        junk = zeros(size(IN));
        [z maxpoint]=max(conv2(r(i,:), w,'same'));
        junk(maxpoint)=1;
        junk = 40*conv2(junk,w,'same');

        all_target_mfr(i,:)=junk;
    end
else
    % Use oriented bar patterns.
    %r=make_patterns(npatterns,length(IN));
    
    r=noise*rand(npatterns,length(IN));
    pat=cos(-1.5*pi:1.5*pi/(Ne/5):1.5*pi);
    pat = pat(pat>0.001);
    pat_width = length(pat);
    spot=1+floor([0 .5 .25 .75]*length(IN));
    gain = signal*[.5 1 .75 .25];
    for i = 1:npatterns
        r(i,spot(i):spot(i)+pat_width-1)=pat*gain(i);
    end

    % DEBUG
    r= rand(npatterns,length(IN));
    r(1,:)=r(1,:)*.5;
    r(2,:)=r(2,:)*1.0;
    r(3,:)=r(3,:)*.75;
    r(4,:)=r(4,:)*.25;

%     r=[.122 .33 .12 .47 1 .42 .03 .43 .31 .31 .45; ...
%     0.2848    0.1522    0.1166    0.2534    0.4529    0.4443    0.1073    0.1756    0.2698    0.3289    0.7500; ...
%     0.1543    0.2126    0.5000    0.1736    0.1401    0.0113    0.0370    0.2063    0.2159    0.2477    0.1655; ...
%     0.1009    0.1317    0.1265    0.0578    0.0548    0.2500    0.1192    0.0408    0.0170    0.0298    0.1115];
%r=[    0.0321    0.4475    0.2013    0.9687    0.1278    0.2822    0.2924    0.1789    0.0484    0.2318    0.7734;...
%    0.7965    0.8375    0.8627    0.0655    0.6800    0.4220    0.9668    0.1903    0.9100    0.5607    0.0139;...
%    0.0341    0.7960    0.2558    0.2295    0.0283    0.6533    0.7794    0.0992    0.1220    0.5733    0.2672;...
%    0.0185    0.9369    0.0597    0.0331    0.8519    0.9216    0.0863    0.9099    0.3537    0.3028    0.9321];

    w=cos(-1.5*pi:1.5*pi/(Ne/4):1.5*pi);
    positive_w=w;
    positive_w(w<0)=0;
    for i = 1:npatterns
        junk = zeros(size(IN));
        [z maxpoint]=max(conv2(r(i,:), w,'same'));
        junk(maxpoint)=1;
        junk = 40*conv2(junk,positive_w,'same');

        all_target_mfr(i,:)=junk;
    end
%     all_target_mfr = [...
%          0.0007    0.0091   14.0048   37.5236   47.7772   37.5828   13.5727    0.0255    0.0067    0.0053    0.0084;...
%          29.3920    7.6146    0.0011    0.0028    0.0073    0.0077    0.0003    0.0069   18.3303   35.9531   42.0335;...
%          18.9614   25.8034   24.2515   11.4287    0.0012         0         0    0.0032    0.0051    0.0224    5.1803;...
%            0.0009    0.0068    0.2328    4.4247    7.4911    8.7554    5.4898    0.3377         0         0    0.0006 ...
%             ];

   
end

for i = 1:npatterns
    patterns{i} = r(i,:);
end



    

%function spikeout = spnet( stop_sec, start_save_sec,  conn_data)
% spikeout will be all spikes between start_savfe_sec and stop_sec.


learning = false;
  

spikeout=[];

% spnet.m: Spiking network with axonal conduction delays and STDP
% Created by Eugene M.Izhikevich.                February 3, 2004
% Modified to allow arbitrary delay distributions.  April 16,2008
% Rewritten by Jeffrey L. McKinstry, November, 2009 to allow anatomy.

D=5;                  % maximal conduction delay 


% neuron parameters from Izhikevich & Edelman (2008) RS and FS cells.
vr=([-60*ones(Ne,1);    -55*ones(Ni,1)]);
vt=([-50*ones(Ne,1);    -40*ones(Ni,1)]);
k=([3*ones(Ne,1);        1*ones(Ni,1)]);
C=([80*ones(Ne,1);      20*ones(Ni,1)]);
a=([0.01*ones(Ne,1);    0.15*ones(Ni,1)]);
b=([5*ones(Ne,1);        8*ones(Ni,1)]);
c=([-60*ones(Ne,1);     -55*ones(Ni,1)]);
d=([400*ones(Ne,1);    200*ones(Ni,1)]);
vp=([50*ones(Ne,1);      25*ones(Ni,1)]);

% changing pyr params to be more like ss in layer IV
group = A;
vr(group)= -60;
vt(group)= -40;
k(group)=   3;
C(group)= 100;
a(group)=0.01;
b(group)=   5; % Assuming it is in tonic mode all the time.
c(group)=   -50;
d(group)=   50;
vp(group)=  35;



group = IN;
vr(group)= -60;
vt(group)= -50;
k(group)=   1.6;
C(group)= 200;
a(group)=0.1;
b(group)=   0; % Assuming it is in tonic mode all the time.
c(group)=   -60;
d(group)=   10;
vp(group)=  40;


% Take special care not to have multiple connections between neurons
s = cell(D,1);
sd = cell(D,1);
syntype = cell(D,1);
for i = 1:N_connections
    for j = 1:D
        s{i,j} = sparse(N,N);
        syntype{i,j} = sparse(N,N);
        sd{i,j} = sparse(N,N);

    end
end


%% connections


syndata = [];


i=1;


% make the initial output the TOTAL output, not max output.
%    std    wtmin,    wt range, radius delaymin, delaymax,lrate, wtmin,wtmax, background ampa rate(Hz), gaba rate(Hz),   nmdagain,
%    gababgain,      f,     d,       u

radius=2;%length(A)/2;
conn_data(i,:)=[0.5*radius 0 80 radius 3 5 .05 0 20   0 0 1 .5 0 0 5 300 .5]; 
if exist('conn_data_in') && ~isempty(conn_data_in)
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(IN,A,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', ...
    conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), ...
    'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9), 'background_ampa_hz', conn(10),...
    'background_gaba_hz', conn(11), 'ampagain', conn(12), 'nmdagain', conn(13),...
    'gabaagain', conn(14), 'gababgain', conn(15),...
    'f', conn(16), 'd', conn(17), 'u', conn(18));


radius = length(A)/2;
not_used = pi;
conn_data(i,:)=[not_used 0  35 radius 1 5 0.0 0 10  0 0 1 0.5 0 0 20 20 .5];
if exist('conn_data_in') && ~isempty(conn_data_in)

    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(A,A,'cos', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', ...
    conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), ...
    'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9), 'background_ampa_hz', conn(10),...
    'background_gaba_hz', conn(11), 'ampagain', conn(12), 'nmdagain', conn(13),...
    'gabaagain', conn(14), 'gababgain', conn(15),...
    'f', conn(16), 'd', conn(17), 'u', conn(18));



radius = length(A)/2;
conn_data(i,:)=[not_used 0 4*80 radius 1 1 0.0 0 20  0 0 0 0 1 .025 5 45 .56]; 
if exist('conn_data_in') && ~isempty(conn_data_in)
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(AI,A,'cossurround', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', ...
    conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), ...
    'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9), 'background_ampa_hz', conn(10),...
    'background_gaba_hz', conn(11), 'ampagain', conn(12), 'nmdagain', conn(13),...
    'gabaagain', conn(14), 'gababgain', conn(15),...
    'f', conn(16), 'd', conn(17), 'u', conn(18));



%radius=length(AI)/2;
radius = 2;
conn_data(i,:)=[0.5*radius 0 20 radius 1 2 .05 0 20  0 0 1 .5 1 0 0 300 .5]; 
if exist('conn_data_in') && ~isempty(conn_data_in)

    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(IN,AI,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', ...
    conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), ...
    'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9), 'background_ampa_hz', conn(10),...
    'background_gaba_hz', conn(11), 'ampagain', conn(12), 'nmdagain', conn(13),...
    'gabaagain', conn(14), 'gababgain', conn(15),...
    'f', conn(16), 'd', conn(17), 'u', conn(18));


radius=2;%0.4;
conn_data(i,:)=[0.5*radius 0 40 radius 1 2 .05 0 20  800 400 1 0.5 0 0 1000 800 .5]; 
%conn_data(i,:)=[0.5*radius 0 80 radius 1 2 .05 0 20]; 
if exist('conn_data_in') && ~isempty(conn_data_in)

    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(A,AI,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', ...
    conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), ...
    'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9), 'background_ampa_hz', conn(10),...
    'background_gaba_hz', conn(11), 'ampagain', conn(12), 'nmdagain', conn(13),...
    'gabaagain', conn(14), 'gababgain', conn(15),...
    'f', conn(16), 'd', conn(17), 'u', conn(18));



radius = length(AI)/2;
conn_data(i,:)=[not_used 0 100 radius 1 1 0.0 0 20  0 0 0 0 1 0 20 700 .2]; 
if exist('conn_data_in') && ~isempty(conn_data_in)
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(AI,AI,'cossurround', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', ...
    conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), ...
    'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9), 'background_ampa_hz', conn(10),...
    'background_gaba_hz', conn(11), 'ampagain', conn(12), 'nmdagain', conn(13),...
    'gabaagain', conn(14), 'gababgain', conn(15),...
    'f', conn(16), 'd', conn(17), 'u', conn(18));



mfr=zeros(1,N);




%s=[6*ones(Ne,M);-5*ones(Ni,M)];         % synaptic weights

if learning
    ltp_decay = exp(-1/20);
    ltpcurve = .1*cumprod([ 1 ltp_decay*ones(1,49)]);
    ltd_decay = exp(-1/20);
    ltdcurve = -.12*cumprod([ 1 ltd_decay*ones(1,49)]);
    stdpcurve = [ fliplr(ltpcurve)  ltdcurve 0];
end  

if draw
    mv = zeros(N,1); % for visualization
    vhist = zeros(N,1000);
    uhist = zeros(N,1000);
    ihist = zeros(N,1000);
    for i = 1:N_connections
        ehist{i} = zeros(N,1000);
    end
    ampahist = zeros(N,1000);
    nmdahist = zeros(N,1000);
    gabaahist = zeros(N,1000);
    gababhist = zeros(N,1000);
    itothist = zeros(N,1000);    
end


% restart for each simulation.
%rand('seed',1)

%% model conductances
tau_ampa  = ones(N,1)*exp(-1/5);
tau_nmda  = ones(N,1)*exp(-1/150);
tau_gabaa = ones(N,1)*exp(-1/6);
tau_gabab = ones(N,1)*exp(-1/150);

tau_gabab(A)=exp(-1/150);

E_AMPA	= 0.0;
E_NMDA	= 0.0;
E_GABAa	= -70;
E_GABAb = -90;

for i = 1:N_connections
    Ie{i}=zeros(N,D);
end

I=zeros(N,N_connections);
ggabaa = I;ggabab = I; gampa = I;gnmda = I;
    



% end model conductances.
%% init stsp variables:

stf = zeros(N,N_connections);
stdep = zeros(N,N_connections);
stsp = zeros(N,N_connections);
for i = 1:N_connections
    stsp(:,i) = syndata(i).U*ones(N,1);
end

% end stsp parameters
%%
v = -65*ones(N,1)+5*rand(N,1);                      % initial values
u = 0.2*ones(N,1)+20*rand(N,1);                             % initial values

MAXSPIKESPERSEC= N/2*1000;
firings=zeros(MAXSPIKESPERSEC,2);

firehist=cell(N,1);


numspikes = 0;

for i=1:N_connections
    snew = [];
    for j=1:D
        snew = [snew s{i,j}];
    end
    %snew= full(snew);
    se{i}=snew;
end


%% simulation
for sec=1:stop_sec                      % simulation of 1 day
  for t=1:1000                          % simulation of 1 sec
    
    %randindex = ceil(N*rand);
    %I(randindex,1)=I(randindex,1)+20;                 % random thalamic input 
%    I(IN,1) = I(IN,1) + 20*patterns{ 1+floor((t-1)/250)};
    
    fired = find(v>=vp);                % indices of fired neurons
    sp_fired = v'>=vp';
    v(fired)=c(fired);  
    u(fired)=u(fired)+d(fired);

    %--------------------------------------
    % Deal with synapses.
    for i = 1:N_connections
    
        sfired = sp_fired.*stsp(:,i)';


        % Short term synaptic plasticity.
        stdep(fired,i)=stdep(fired,i)+stsp(fired,i);
        stf(fired,i)=stf(fired,i) + syndata(i).U *(1-syndata(i).U-stf(fired,i));

        stf(:,i) = stf(:,i) * syndata(i).tau_F1;
        stdep(:,i) = stdep(:,i)*syndata(i).tau_D1;
        stsp(:,i) = (1-stdep(:,i)).*(syndata(i).U+stf(:,i)); % By default, keep the value of the last spike?

        
        Ie{i} = Ie{i} + reshape( sfired*se{i} ,N,D);

        
        ggabaa(:,i) = ggabaa(:,i).*tau_gabaa + syndata(i).gain_gabaa.*Ie{i}(:,1) ; 
        ggabab(:,i) = ggabab(:,i).*tau_gabab + syndata(i).gain_gabab.*Ie{i}(:,1); % This one should build slowly, but continue
            % to build with only a brief pulse.

        gampa(:,i) = gampa(:,i).*tau_ampa + syndata(i).gain_ampa.*Ie{i}(:,1);
        gnmda(:,i) = gnmda(:,i).*tau_nmda + syndata(i).gain_nmda.*Ie{i}(:,1); % should delay an extra ms for nmda.

        % shift the buffer now.
    
        if draw,
            ehist{i}(:,t)= Ie{i}(:,1); 
        end
        Ie{i} = [Ie{i}(:,2:D) zeros(N,1)];  % rotate.     

               
        
    end
    
            % Combine the conductances from all sources.
    Igabaa = sum(ggabaa,2);
    Igabab = sum(ggabab,2);
    Iampa  = sum(gampa,2);
    Inmda  = sum(gnmda,2);

    
    
    % Done dealing with synapses.
    %--------------------------------------
   
    xx = (-80-v)/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
%    NMDAgate = 1.0;
    
    % DEBUG
%     w=cos(-1.5*pi:1.5*pi/(Ne/4):1.5*pi)';
%     w1=-w;w1(w1<0)=0;
%     w(w<0)=0;
%     w=[zeros(size(IN))' ;w(1:Ne/2); zeros(Ni,1)];
%     w1=[zeros(size(IN))' ;w1(1:Ne/2); zeros(Ni,1)];
%     gampa_rand=w*5.*abs(randn(size(w)));
%     ggabaa_rand=w1*20.*abs(randn(size(w)));
    
    %gampa_rand  = abs(randn(N,1)).*gampa_background;
    %ggabaa_rand = abs(randn(N,1)).*ggabaa_background;

    g= ( Iampa+ NMDAgate.*Inmda        +Igabaa+ Igabab);
	E= ( (Iampa)*E_AMPA             +NMDAgate.*Inmda*E_NMDA +(Igabaa)*E_GABAa +Igabab*E_GABAb);
%   g= ( Iampa+ gampa_rand  +NMDAgate.*Inmda        +Igabaa+ ggabaa_rand
%    +Igabab);%
%	E= ( (Iampa+gampa_rand)*E_AMPA             +NMDAgate.*Inmda*E_NMDA +(Igabaa+ggabaa_rand)*E_GABAa +Igabab*E_GABAb);
    Isyn = v.*g-E;
    
    % gap junctions
%     Isyn(AI(2:end-1)) = Isyn(AI(2:end-1)) + (-v(AI(1:end-2))+v(AI(2:end-1)))*6 ...
%                         + (-v(AI(3:end))+v(AI(2:end-1)))*6;
%     Isyn(A(2:end-1)) = Isyn(A(2:end-1)) + (-v(A(1:end-2))+v(A(2:end-1)))*6 ...
%                         + (-v(A(3:end))+v(A(2:end-1)))*6;
    
    %E(IN)=E(IN)+ 40*patterns{ 1+floor((t-1)/isi_ms)}';
	%Isyn(A) = Isyn(A) + -2000*patterns{ 1+floor((t-1)/isi_ms)}';
	Isyn(IN) = -1000*patterns{ 1+floor((t-1)/isi_ms)}';
	%Isyn(IN) = -100*ceil(t/100);
    
    
    v=v+0.5*(k.*(v-vr).*(v-vt)-u-Isyn)./C;    % for numerical 
    %v(v>104)=104;   % Nothing higher than Eca2+
    v(v>50)=50;   % Nothing higher than typical action potential
    v(v<-90)=-90;  % cheating on numerical integration. Nothing lower than EK

    u=u+0.5*a.*(b.*(v-vr)-u);                   % step is 0.5 ms
 
    v=v+0.5*(k.*(v-vr).*(v-vt)-u-Isyn)./C;    % for numerical 
    %v(v>104)=104;   % Nothing higher than Eca2+
    v(v>50)=50;   % Nothing higher than Eca2+
    v(v<-90)=-90;  % cheating on numerical integration. Nothing lower than EK

    u=u+0.5*a.*(b.*(v-vr)-u);                   % step is 0.5 ms
    
    %v(AI)=conv2(v(AI),[.05 .1 .7 .1 .05]','same');
    %v(A)=conv2(v(A),[.1 .8 .1]','same');


    

    
    if draw, 
        itothist(:,t) = Isyn;
        ampahist(:,t)= Iampa; nmdahist(:,t)= Inmda; 
        gabaahist(:,t)=Igabaa; gababhist(:,t)=Igabab;
    end
    
    
    % save the spike data.
    firings(numspikes+1:numspikes+length(fired),:)=[t*ones(length(fired),1),fired];
    numspikes = numspikes + length(fired);
    
    if learning
        for i =1:length(fired)
            firehist{fired(i)}(end+1)=t;
        end
    end
    %imagesc(I);
  
    if animate && mod(t,5)==0
        figure(3)
        for i = 1:N_connections
          subplot(N_connections,1,i)
          plot([gampa(syndata(i).post,i) gnmda(syndata(i).post,i) ggabaa(syndata(i).post,i) ggabab(syndata(i).post,i)]);
          if i==1
              legend('ampa','nmda','gabaa','gabab');
          end
          axis([1 length(A) 0 40]);
          title(conn_name{i});
        end
        text(length(A)/2,-10,num2str(t));
        drawnow
    end
    
    
    
    if draw
        vhist(:,t)=v;
        uhist(:,t)=u;
        
        tau = exp(-1/200);
        mfr(fired)=mfr(fired) + (1-tau)*1000;
        mfr=mfr*tau;

        mv = tau*mv + (1-tau)*v;
        if mod(t,250)==0
            figure(3);
            bar(mfr)
            drawnow;
            
            %figure(4)
            %plot(mv);
            %drawnow;
        end
    end
    
  end;
  
  
  disp(['t=' num2str(sec)])
  
  if draw
      
      figure(1)
      plot(firings(:,1),firings(:,2),'.');
      axis([0 1000 0 N]); drawnow;

      %figure(2)
      %hist(s{1}(find(s{1})),50);

      figure(2)
      for i = 1:N_connections
          subplot(N_connections,1,i)
          imagesc(ehist{i}); colorbar;
          title(conn_name{i});
      end
      
      
      figure(4)
      subplot(7,1,1)
      imagesc(vhist,[-70 -60]); colorbar;
      title('vhist')
      subplot(7,1,2)
      imagesc(ihist); colorbar;
      title('inhib input history');
      subplot(7,1,4)
      imagesc(ampahist); colorbar;
      title('ampa input history');
      subplot(7,1,5)
      imagesc(nmdahist); colorbar;
      title('nmda input history');
      subplot(7,1,6)
      imagesc(gabaahist); colorbar;
      title('gabaa input history');
      subplot(7,1,7)
      imagesc(gababhist); colorbar;
      title('gabab input history');
      
      figure(5)
      subplot(3,1,1)
      j=1;
      for i = 1:isi_ms:999
          sum_temp(:,j) = mean(itothist(:,i:i+99),2)  ;
          %sum_tempe(:,j) = mean(ehist(:,i:i+99),2) ;
          sum_tempv(:,j) = mean(vhist(:,i:i+99),2) ;
          j=j+1;
      end
      plot(sum_temp);
      size(sum_temp)
      title('Isyn')
      
      legend('1','2','3','4','5','6','7','8');
      subplot(3,1,2)
      %plot(sum_tempe);
      title('ehist');

      subplot(3,1,3)
      plot(sum_tempv);
      title('voltage hist');
      
  end
  % do some learning based on spike times.
  % sd in Eugene's spnet did not get used until after 1 second, so this is
  % equivalent.  Just makes for a cleaner implementation.
  
  
  if learning
      for delay = 1:D
        % get all the nonzero elements from synaptic matrix.
        [pre_i post_i] = find( s{delay} );

        % Adjust these weights.
        for j = 1:length(pre_i)
            type = syntype{delay}(pre_i(j),post_i(j));
            w=s{delay}(pre_i(j),post_i(j));

            %if type == 6, w,  sd{delay}(pre_i(j),post_i(j)),end

            if syndata(type).max < .00001
                scale = 1;
            else
                scale = w/ syndata(type).max;
            end
            sd{delay}(pre_i(j),post_i(j)) = sd{delay}(pre_i(j),post_i(j)) + stdp(firehist{pre_i(j)},firehist{post_i(j)}, delay, stdpcurve,syndata(type).lrate, scale);

            %if type == 6, w,  sd{delay}(pre_i(j),post_i(j)),end
            s{delay}(pre_i(j),post_i(j)) = max( syndata(type).min ,min(syndata(type).max, syndata(type).lrate*.01 + w + sd{delay}(pre_i(j),post_i(j)))); 
        end

        sd{delay}=0.9*sd{delay};

      end
      % keep spikes near the seconds boundary for future STDP calculations.
      for i = 1:N
        firehist{i}=firehist{i}( find(firehist{i}> 950) )-1000;
      end
      
  end
  
  if sec >= start_save_sec
      spikeout = [spikeout; [firings(1:numspikes,1)+1000*(sec-1) firings(1:numspikes,2)]];
  end
  
  numspikes = 0;
  
  
end;


if draw
    figure(7)
    %close all
    subplot(2,1,1)
    m=A(20);
    plot([vhist(m,:)' ampahist(m,:)' gabaahist(m,:)' itothist(m,:)' ])
    legend('vhist(m,:)', 'ampahist(m,:)', 'gabaahist(m,:)', 'itothist(m,:)' )
    
    subplot(2,1,2)
    plot(spikeout(:,1),spikeout(:,2),'.');
    
end


% Calculate cost measure

% build the desired mfr profile.


% all_target_mfr = zeros(npatterns,length(A));
% for pat = 1:npatterns
% 
%     std=length(A)/10;
%     radius=floor(length(A)/2);
%     target_mfr=(gaussian_1d(std,-radius:radius));
%     target_mfr=target_mfr*(1/max(target_mfr)).*30; % winner fires at 30 hz.
%     [m max_index]=max( patterns{pat} );  % shift s.t. winner is in the middle.
%     target_mfr=target_mfr(1:length(A));
%     % assumes length(A) = length(r);
%     mid = ceil(length(A)/2);
%     shift = abs(max_index - mid)
%     if max_index > mid
%         % shift right
%         target_mfr=[ target_mfr(end-shift+1:end) target_mfr(1:end-shift)];
%     else
%         % mid >= max_index
%         target_mfr=[ target_mfr(shift+1:end) target_mfr(1:shift)];    
%     end
% 
%     all_target_mfr(pat,:)=target_mfr;
% end
    
% bin the spike data by pattern.
if ~isempty(spikeout)    
    S=sparse( spikeout(:,2), spikeout(:,1), 1, N , stop_sec*1000 );

    % Gather mean rates for stats.

    t=start_save_sec*1000 + 1;
    count = zeros(1,npatterns);
    mfr = zeros(npatterns,N);
    for i = start_save_sec*1000:epoch_ms:stop_sec*1000-1
        for p = 1:npatterns
            mfr(p, : ) = mfr(p, :) + (1000/isi_ms)*full(sum( S( :, t+20:(t+isi_ms-1)), 2))';
            count(p)=count(p)+1;
            t=t+isi_ms;
        end
    end
    mfr=mfr./(repmat( count', 1, size(mfr,2) ) + .000001);


    % Look at sparsness measures
    %     
    %     j= 1;
    %     for i = max(IN)+1:Ne
    %         lifetime_sparseness(j) = sparseness( mfr(:,i), 40 ); % 40 Hz is max
    %         j= j + 1;
    %     end
    %     disp(['Mean lifetime sparseness: ' num2str( mean(lifetime_sparseness))]);
    %     
    %     for i = 1:patterns
    %         population_sparseness(i) = sparseness(mfr(i,max(IN)+1:Ne), 40);
    %     end
    %     disp(['Mean population sparseness: ' num2str( mean(population_sparseness))]);

        % save the best
    %    d= sqrt( (mean(population_sparseness)-0.5)^2 )
    %d= sqrt( (mean(population_sparseness)-0.5)^2 + (mean(lifetime_sparseness)-0.5)^2)

    cost_measure = sqrt( sum(sum((mfr(:,max(IN)+1:Ne) - all_target_mfr ).^2 )))

else
    cost_measure = 1e10;
end


if draw

    figure(8)
    plot(mfr(:,A)','--')
    hold on;plot(all_target_mfr')
    
    figure(9)
    subplot(2,1,1)
    plot(ampahist(A(20),:))
    hold on
    plot(nmdahist(A(20),:),'g')
    plot(vhist(20,:),'r')
    hold off
    subplot(2,1,2)
    plot(vhist(A(20),:))
    
end



%
% possible arguments:
%  arbor_type == 'gauss'
%  'std', 'minw', 'maxrangew', 'dist', 'delayrange', 'lrate';
%
function connect( from_indices, to_indices, arbor_type, varargin)

global s sd syntype N D syndata


L = length( syndata);

syndata(L+1).post=to_indices;
syndata(L+1).pre =from_indices;

for k= 1 :2: size(varargin,2) 

        switch lower( varargin{k} )
          case {'std'}
            std = varargin{k+1};
          case 'initoffsetw'
            minw = varargin{k+1};
          case 'initrangew'
            maxw = varargin{k+1};
          case 'radius'
            radius = round(varargin{k+1});
          case 'delaymin'
            delaymin = round(varargin{k+1});
          case 'delayrange'
            delayrange = round(varargin{k+1});
            if delaymin + delayrange-1 > D
                disp('delay > max delay in connect. Using max.');
                delayrange
                delaymin
                delayrange = D-delaymin + 1;
            end
          case 'lrate'
            syndata(L+1).lrate = varargin{k+1};
          case 'minw'
            syndata(L+1).min = varargin{k+1};
          case 'maxrangew'
            syndata(L+1).max = syndata(L+1).min + varargin{k+1};
          case 'background_ampa_hz'
            syndata(L+1).background_ampa_hz = varargin{k+1};
          case 'background_gaba_hz'
            syndata(L+1).background_gaba_hz = varargin{k+1};
          case 'ampagain'
            syndata(L+1).gain_ampa = varargin{k+1};
          case 'nmdagain'
            syndata(L+1).gain_nmda = varargin{k+1};
          case 'gabaagain'
            syndata(L+1).gain_gabaa = varargin{k+1};
          case 'gababgain'
            syndata(L+1).gain_gabab = varargin{k+1};
          case 'f'
            syndata(L+1).tau_F1=exp(-1/varargin{k+1});
          case 'd'
            syndata(L+1).tau_D1=exp(-1/varargin{k+1});
          case 'u'
            syndata(L+1).U = varargin{k+1};



          otherwise
            disp(['Unknown param.' varargin{k}])
        end

end

% Assume wraparound.

% clip radius to not touch target neurons more than once.
if radius > length(to_indices)/2-.5, radius = floor(length(to_indices)/2-.5);end


p=[to_indices to_indices to_indices]; %randperm(N);
% create the list of indices for each neuron; project topographically.
step = length(to_indices)/length(from_indices);
rel_index = .5;

if strcmp(arbor_type, 'gauss')
    w = std*sqrt(2*pi)* gaussian_1d(std, -round(radius):round(radius) );
    w = minw+maxw *(1/sum(w))*w;
elseif strcmp(arbor_type, 'surround')
    % surround
    x=-round(radius):round(radius);
    w=std*sqrt(2*pi)* abs(x).*gaussian_1d(std, x );
    w = minw+maxw*(1/sum(w))*w;
elseif strcmp(arbor_type, 'cos')
    w=cos(-1.5*pi:1.5*pi/radius:1.5*pi);
    w(w<0)=0;
    w=minw+maxw*(1/sum(w))*w;  % make them sum to maxw.
elseif strcmp(arbor_type, 'cossurround')
    w=cos(-1.5*pi:1.5*pi/radius:1.5*pi);
    w(w>0)=0;
    w=-w; % make these weights positive for the conductance-based model.
    w=minw+maxw*(1/sum(w))*w;  % make them sum to maxw.
    
end
for i=from_indices
    my_center = round(rel_index*step);
    rel_index=rel_index + 1;
    myoffset = length(to_indices) + my_center;

    target_indices = p( myoffset + (-round(radius):round(radius)) );
    
    delays = delaymin + floor( delayrange*(abs(-round(radius):round(radius))/(round(radius)+.001)));

    for j = 1:length(target_indices)
        s{L+1,delays(j)}(i,target_indices(j)) = w(j); % + noise.
        syntype{L+1,delays(j)}(i,target_indices(j)) = length( syndata );
        sd{L+1,delays(j)}(i,target_indices(j)) = 0.0000001;         
    end

    
end;


function xa = gaussian_1d(sigma_p,x);

x0 = 0;


xa = 1/(sigma_p * sqrt(2*pi)) * exp ((-(-x-x0).^2) ./ (2*sigma_p^2));




%%
% make some oriented gabor patches.
function P = make_patterns(samples,n);


pixels = floor(floor(sqrt(n/2)));


P = zeros( samples,n);
for i = 1:samples
    cycles = 2;
    xphase = 0; %(-1+2*rand)*cycles*pi;
    yphase = 0; %(-1+2*rand)*cycles*pi;
    angle = i*pi/samples;
    contrast = .25+.75*rand;
    sd_x = 1;
    sd_y=4;

    g=contrast * gaborfilter(pixels,sd_x*pi,sd_y*pi,angle,cycles,xphase,yphase);

%     if mod(i,1)==0
%         imagesc(g,[-1 1])
%         colormap(gray);
%         drawnow;
%     end

    on=g(:);
    off=-on;
    
    on(on<0) = 0;
    off(off<0)=0;
    temp=[on;off];    
    P(i,1:length(temp))= temp';

    %pause(.5);    
    
end

